#undef DEBUG
module cldwat
!----------------------------------------------------------------------- 
! 
! Purpose: Prognostic cloud water data and methods.
! 
! Public interfaces:
!
! inimc -- Initialize constants
! pcond -- Calculate prognostic condensate
!
! cldwat_fice   -- calculate fraction of condensate in ice phase (radiation partitioning)
!
! Author: P. Rasch, with Modifications by Minghua Zhang
! January 2010, modified by J. Kay to add precip fluxes for COSP simulator
! 
!-----------------------------------------------------------------------
   use shr_kind_mod,  only: r8 => shr_kind_r8
   use spmd_utils,    only: masterproc
   use ppgrid,        only: pcols, pver, pverp
   use wv_saturation, only: estblf, hlatv, tmin, hlatf, rgasv, pcf, &
                            cp, epsqs, ttrice
   use abortutils,    only: endrun
   use cam_logfile,   only: iulog

   implicit none

!-----------------------------------------------------------------------
! PUBLIC: Make default data and interfaces private
!-----------------------------------------------------------------------
   private
   save
   public inimc, pcond          ! Public interfaces
   public cldwat_fice          ! Public interfaces
   public cldwat_readnl
   integer, public::  ktop      ! Level above 10 hPa

   real(r8),public ::  icritc               ! threshold for autoconversion of cold ice
   real(r8),public ::  icritw               ! threshold for autoconversion of warm ice
!!$   real(r8),public,parameter::  conke  = 1.e-6    ! tunable constant for evaporation of precip
!!$   real(r8),public,parameter::  conke  =  2.e-6    ! tunable constant for evaporation of precip
   real(r8),public ::  conke                ! tunable constant for evaporation of precip
   real(r8),public ::  r3lcrit              ! critical radius where liq conversion begins

!-----------------------------------------------------------------------
! PRIVATE: Everything else is private to this module
!-----------------------------------------------------------------------
   real(r8), private:: tmax_fice! max temperature for cloud ice formation
   real(r8), private:: tmin_fice! min temperature for cloud ice formation
   real(r8), private:: tmax_fsnow ! max temperature for transition to convective snow
   real(r8), private:: tmin_fsnow ! min temperature for transition to convective snow
   real(r8), private:: rhonot   ! air density at surface
   real(r8), private:: t0       ! Freezing temperature
   real(r8), private:: cldmin   ! assumed minimum cloud amount
   real(r8), private:: small    ! small number compared to unity
   real(r8), private:: c        ! constant for graupel like snow cm**(1-d)/s
   real(r8), private:: d        ! constant for graupel like snow
   real(r8), private:: esi      ! collection efficient for ice by snow
   real(r8), private:: esw      ! collection efficient for water by snow
   real(r8), private:: nos      ! particles snow / cm**4
   real(r8), private:: pi       ! Mathematical constant
   real(r8), private:: gravit   ! Gravitational acceleration at surface
   real(r8), private:: rh2o
   real(r8), private:: prhonos
   real(r8), private:: thrpd    ! numerical three added to d
   real(r8), private:: gam3pd   ! gamma function on (3+d)
   real(r8), private:: gam4pd   ! gamma function on (4+d)
   real(r8), private:: rhoi     ! ice density
   real(r8), private:: rhos     ! snow density
   real(r8), private:: rhow     ! water density
   real(r8), private:: mcon01   ! constants used in cloud microphysics
   real(r8), private:: mcon02   ! constants used in cloud microphysics
   real(r8), private:: mcon03   ! constants used in cloud microphysics
   real(r8), private:: mcon04   ! constants used in cloud microphysics
   real(r8), private:: mcon05   ! constants used in cloud microphysics
   real(r8), private:: mcon06   ! constants used in cloud microphysics
   real(r8), private:: mcon07   ! constants used in cloud microphysics
   real(r8), private:: mcon08   ! constants used in cloud microphysics

   integer, private ::  k1mb    ! index of the eta level near 1 mb

! Parameters used in findmcnew
   real(r8) :: capnsi               ! sea ice cloud particles / cm3
   real(r8) :: capnc                ! cold and oceanic cloud particles / cm3
   real(r8) :: capnw                ! warm continental cloud particles / cm3
   real(r8) :: kconst               ! const for terminal velocity (stokes regime)
   real(r8) :: effc                 ! collection efficiency
   real(r8) :: alpha                ! ratio of 3rd moment radius to 2nd
   real(r8) :: capc                 ! constant for autoconversion
   real(r8) :: convfw               ! constant used for fall velocity calculation
   real(r8) :: cracw                ! constant used for rain accreting water
   real(r8) :: critpr               ! critical precip rate collection efficiency changes
   real(r8) :: ciautb               ! coefficient of autoconversion of ice (1/s)

#ifdef DEBUG
   integer, private,parameter ::  nlook = 1  ! Number of points to examine
   integer, private ::  ilook(nlook)         ! Longitude index to examine
   integer, private ::  latlook(nlook)       ! Latitude index to examine
   integer, private ::  lchnklook(nlook)     ! Chunk index to examine
   integer, private ::  icollook(nlook)      ! Column index to examine
#endif
  ! Private data
  real(r8), parameter :: unset_r8 = huge(1.0_r8)

contains
!===============================================================================
  subroutine cldwat_readnl(nlfile)

   use namelist_utils,  only: find_group_name
   use units,           only: getunit, freeunit
   use mpishorthand

   character(len=*), intent(in) :: nlfile  ! filepath for file containing namelist input




   ! Namelist variables
   real(r8) :: cldwat_icritw  = unset_r8    !   icritw  = threshold for autoconversion of warm ice  
   real(r8) :: cldwat_icritc  = unset_r8    !   icritc  = threshold for autoconversion of cold ice  
   real(r8) :: cldwat_conke   = unset_r8    !   conke   = tunable constant for evaporation of precip
   real(r8) :: cldwat_r3lcrit = unset_r8    !   r3lcrit = critical radius where liq conversion begins

   ! Local variables
   integer :: unitn, ierr
   character(len=*), parameter :: subname = 'cldwat_readnl'

   namelist /cldwat_nl/ cldwat_icritw, cldwat_icritc, cldwat_conke, cldwat_r3lcrit

   !-----------------------------------------------------------------------------

   if (masterproc) then
      unitn = getunit()
      open( unitn, file=trim(nlfile), status='old' )
      call find_group_name(unitn, 'cldwat_nl', status=ierr)
      if (ierr == 0) then
         read(unitn, cldwat_nl, iostat=ierr)
         if (ierr /= 0) then
            call endrun(subname // ':: ERROR reading namelist')
         end if
      end if
      close(unitn)
      call freeunit(unitn)

      ! set local variables
      icritw  = cldwat_icritw 
      icritc  = cldwat_icritc
      conke   = cldwat_conke
      r3lcrit = cldwat_r3lcrit

   end if



#ifdef SPMD
   ! Broadcast namelist variables
   call mpibcast(icritw,            1, mpir8,  0, mpicom)
   call mpibcast(icritc,            1, mpir8,  0, mpicom)
   call mpibcast(conke,             1, mpir8,  0, mpicom)
   call mpibcast(r3lcrit,           1, mpir8,  0, mpicom)
#endif

end subroutine cldwat_readnl

!================================================================================================
  subroutine cldwat_fice(ncol, t, fice, fsnow)
!
! Compute the fraction of the total cloud water which is in ice phase.
! The fraction depends on temperature only. 
! This is the form that was used for radiation, the code came from cldefr originally
! 
! Author: B. A. Boville Sept 10, 2002
!  modified: PJR 3/13/03 (added fsnow to ascribe snow production for convection )
!-----------------------------------------------------------------------
    implicit none

! Arguments
    integer,  intent(in)  :: ncol                 ! number of active columns
    real(r8), intent(in)  :: t(pcols,pver)        ! temperature

    real(r8), intent(out) :: fice(pcols,pver)     ! Fractional ice content within cloud
    real(r8), intent(out) :: fsnow(pcols,pver)    ! Fractional snow content for convection

! Local variables
    integer :: i,k                                   ! loop indexes

!-----------------------------------------------------------------------

! Define fractional amount of cloud that is ice
    do k=1,pver
       do i=1,ncol

! If warmer than tmax then water phase
          if (t(i,k) > tmax_fice) then
             fice(i,k) = 0.0_r8

! If colder than tmin then ice phase
          else if (t(i,k) < tmin_fice) then
             fice(i,k) = 1.0_r8

! Otherwise mixed phase, with ice fraction decreasing linearly from tmin to tmax
          else 
             fice(i,k) =(tmax_fice - t(i,k)) / (tmax_fice - tmin_fice)
          end if

! snow fraction partitioning

! If warmer than tmax then water phase
          if (t(i,k) > tmax_fsnow) then
             fsnow(i,k) = 0.0_r8

! If colder than tmin then ice phase
          else if (t(i,k) < tmin_fsnow) then
             fsnow(i,k) = 1.0_r8

! Otherwise mixed phase, with ice fraction decreasing linearly from tmin to tmax
          else 
             fsnow(i,k) =(tmax_fsnow - t(i,k)) / (tmax_fsnow - tmin_fsnow)
          end if

       end do
    end do

    return
  end subroutine cldwat_fice

subroutine inimc( tmeltx, rhonotx, gravitx, rh2ox)
!----------------------------------------------------------------------- 
! 
! Purpose: 
! initialize constants for the prognostic condensate
! 
! Author: P. Rasch, April 1997
! 
!-----------------------------------------------------------------------
   use pmgrid,       only: plev, plevp
   use dycore,       only: dycore_is, get_resolution
   use hycoef,       only: hypm
   use phys_control, only: phys_getopts

   integer k
   real(r8), intent(in) :: tmeltx
   real(r8), intent(in) :: rhonotx
   real(r8), intent(in) :: gravitx
   real(r8), intent(in) :: rh2ox

#ifdef UNICOSMP
   real(r8) signgam              ! variable required by cray gamma function
   external gamma
#endif
   character(len=16)    :: microp_scheme 

! Get microphysics option

   call phys_getopts( microp_scheme_out = microp_scheme )

! Set following for all physics packages

   tmax_fice = tmeltx    - 10._r8
!! tmax_fice = tmeltx
!! tmin_fice = tmax_fice - 20.
   tmin_fice = tmax_fice - 30._r8
   tmax_fsnow = tmeltx
   tmin_fsnow = tmeltx   - 5._r8

! Set remaining for RK microphysics

   if( microp_scheme .eq. 'RK' ) then
      rhonot = rhonotx             ! air density at surface (gm/cm3)
      gravit = gravitx
      rh2o   = rh2ox
      rhos = .1_r8                 ! assumed snow density (gm/cm3)
      rhow = 1._r8                 ! water density
      rhoi = 1._r8                 ! ice density
      esi = 1.0_r8                 ! collection efficient for ice by snow
      esw = 0.1_r8                 ! collection efficient for water by snow
      t0 = tmeltx                  ! approximate freezing temp
      cldmin = 0.02_r8             ! assumed minimum cloud amount
      small = 1.e-22_r8            ! a small number compared to unity
      c = 152.93_r8                ! constant for graupel like snow cm**(1-d)/s
      d = 0.25_r8                  ! constant for graupel like snow
      nos = 3.e-2_r8               ! particles snow / cm**4
      pi = 4._r8*atan(1.0_r8)
      prhonos = pi*rhos*nos
      thrpd = 3._r8 + d
      if (d==0.25_r8) then
         gam3pd = 2.549256966718531_r8 ! only right for d = 0.25
         gam4pd = 8.285085141835282_r8
      else
#ifdef UNICOSMP
         call gamma(3._r8+d, signgam, gam3pd)
         gam3pd = sign(exp(gam3pd),signgam)
         call gamma(4._r8+d, signgam, gam4pd)
         gam4pd = sign(exp(gam4pd),signgam)
         write(iulog,*) ' d, gamma(3+d), gamma(4+d) =', gam3pd, gam4pd
#else
         write(iulog,*) ' can only use d ne 0.25 on a cray '
         stop
#endif
      endif
      mcon01 = pi*nos*c*gam3pd/4._r8
      mcon02 = 1._r8/(c*gam4pd*sqrt(rhonot)/(6*prhonos**(d/4._r8)))
      mcon03 = -(0.5_r8+d/4._r8)
      mcon04 = 4._r8/(4._r8+d)
      mcon05 = (3+d)/(4+d)
      mcon06 = (3+d)/4._r8
      mcon07 = mcon01*sqrt(rhonot)*mcon02**mcon05/prhonos**mcon06
      mcon08 = -0.5_r8/(4._r8+d)

!  find the level about 1mb, we wont do the microphysics above this level
      k1mb = 1
      do k=1,pver-1
         if (hypm(k) < 1.e2_r8 .and. hypm(k+1) >= 1.e2_r8) then
            if (1.e2_r8-hypm(k) < hypm(k+1)-1.e2_r8) then
               k1mb = k
            else
               k1mb = k + 1
            end if
            goto 20
         end if
      end do
      if (masterproc) then
         write(iulog,*)'inimc: model levels bracketing 1 mb not found'
      end if
!     call endrun
      k1mb = 1
20    if( masterproc ) write(iulog,*)'inimc: model level nearest 1 mb is',k1mb,'which is',hypm(k1mb),'pascals'

      if( masterproc ) write(iulog,*) 'cloud water initialization by inimc complete '

! Initialize parameters used by findmcnew
      capnw = 400._r8              ! warm continental cloud particles / cm3
      capnc = 150._r8              ! cold and oceanic cloud particles / cm3
!     capnsi = 40._r8              ! sea ice cloud particles density  / cm3
      capnsi = 75._r8              ! sea ice cloud particles density  / cm3

      kconst = 1.18e6_r8           ! const for terminal velocity

!     effc = 1._r8                 ! autoconv collection efficiency following boucher 96
!     effc = .55*0.05_r8           ! autoconv collection efficiency following baker 93
      effc = 0.55_r8               ! autoconv collection efficiency following tripoli and cotton
!     effc = 0._r8    ! turn off warm-cloud autoconv
      alpha = 1.1_r8**4
      capc = pi**(-.333_r8)*kconst*effc *(0.75_r8)**(1.333_r8)*alpha  ! constant for autoconversion

! critical precip rate at which we assume the collector drops can change the
! drop size enough to enhance the auto-conversion process (mm/day)
      critpr = 0.5_r8
      convfw = 1.94_r8*2.13_r8*sqrt(rhow*1000._r8*9.81_r8*2.7e-4_r8)

! liquid microphysics
!     cracw = 6_r8                 ! beheng
      cracw = .884_r8*sqrt(9.81_r8/(rhow*1000._r8*2.7e-4_r8)) ! tripoli and cotton

! ice microphysics
      ciautb = 5.e-4_r8

      if ( masterproc ) then
         write(iulog,*)'tuning parameters cldwat: icritw',icritw,'icritc',icritc,'conke',conke,'r3lcrit',r3lcrit
         write(iulog,*)'tuning parameters cldwat: capnw',capnw,'capnc',capnc,'capnsi',capnsi,'kconst',kconst
         write(iulog,*)'tuning parameters cldwat: effc',effc,'alpha',alpha,'capc',capc
         write(iulog,*)'tuning parameters cldwat: critpr',critpr,'convfw',convfw,'cracw',cracw,'ciautb',ciautb
      endif
   endif

   return
end subroutine inimc

subroutine pcond (lchnk   ,ncol    , &
                  tn      ,ttend   ,qn      ,qtend   ,omega   , &
                  cwat    ,p       ,pdel    ,cldn    ,fice    , fsnow, &
                  cme     ,prodprec,prodsnow,evapprec,evapsnow,evapheat, prfzheat, &     
                  meltheat,precip  ,snowab  ,deltat  ,fwaut   , &
                  fsaut   ,fracw   ,fsacw   ,fsaci   ,lctend  , &
                  rhdfda  ,rhu00   ,seaicef, zi      ,ice2pr, liq2pr, &
                  liq2snow, snowh, rkflxprc, rkflxsnw, pracwo, psacwo, psacio)

!----------------------------------------------------------------------- 
! 
! Purpose: 
! The public interface to the cloud water parameterization
! returns tendencies to water vapor, temperature and cloud water variables
! 
! For basic method 
!  See: Rasch, P. J, and J. E. Kristjansson, A Comparison of the CCM3
!  model climate using diagnosed and 
!  predicted condensate parameterizations, 1998, J. Clim., 11,
!  pp1587---1614.
! 
! For important modifications to improve the method of determining
! condensation/evaporation see Zhang et al (2001, in preparation)
!
! Authors: M. Zhang, W. Lin, P. Rasch and J.E. Kristjansson
!          B. A. Boville (latent heat of fusion)
!-----------------------------------------------------------------------
   use wv_saturation, only: vqsatd
   use cam_control_mod, only: nlvdry
!
!---------------------------------------------------------------------
!
! Input Arguments
!
   integer, intent(in) :: lchnk                 ! chunk identifier
   integer, intent(in) :: ncol                  ! number of atmospheric columns

   real(r8), intent(in) :: fice(pcols,pver)     ! fraction of cwat that is ice
   real(r8), intent(in) :: fsnow(pcols,pver)    ! fraction of rain that freezes to snow
   real(r8), intent(in) :: cldn(pcols,pver)     ! new value of cloud fraction    (fraction)
   real(r8), intent(in) :: cwat(pcols,pver)     ! cloud water (kg/kg)
   real(r8), intent(in) :: omega(pcols,pver)    ! vert pressure vel (Pa/s)
   real(r8), intent(in) :: p(pcols,pver)        ! pressure          (K)
   real(r8), intent(in) :: pdel(pcols,pver)     ! pressure thickness (Pa)
   real(r8), intent(in) :: qn(pcols,pver)       ! new water vapor    (kg/kg)
   real(r8), intent(in) :: qtend(pcols,pver)    ! mixing ratio tend  (kg/kg/s)
   real(r8), intent(in) :: tn(pcols,pver)       ! new temperature    (K)
   real(r8), intent(in) :: ttend(pcols,pver)    ! temp tendencies    (K/s)
   real(r8), intent(in) :: deltat               ! time step to advance solution over
   real(r8), intent(in) :: lctend(pcols,pver)   ! cloud liquid water tendencies   ====wlin
   real(r8), intent(in) :: rhdfda(pcols,pver)   ! dG(a)/da, rh=G(a), when rh>u00  ====wlin
   real(r8), intent(in) :: rhu00 (pcols,pver)   ! Rhlim for cloud                 ====wlin
   real(r8), intent(in) :: seaicef(pcols)       ! sea ice fraction  (fraction)
   real(r8), intent(in) :: zi(pcols,pverp)      ! layer interfaces (m)
    real(r8), intent(in) :: snowh(pcols)         ! Snow depth over land, water equivalent (m)
!
! Output Arguments
!
   real(r8), intent(out) :: cme     (pcols,pver) ! rate of cond-evap of condensate (1/s)
   real(r8), intent(out) :: prodprec(pcols,pver) ! rate of conversion of condensate to precip (1/s)
   real(r8), intent(out) :: evapprec(pcols,pver) ! rate of evaporation of falling precip (1/s)
   real(r8), intent(out) :: evapsnow(pcols,pver) ! rate of evaporation of falling snow (1/s)
   real(r8), intent(out) :: evapheat(pcols,pver) ! heating rate due to evaporation of precip (W/kg)
   real(r8), intent(out) :: prfzheat(pcols,pver) ! heating rate due to freezing of precip (W/kg)
   real(r8), intent(out) :: meltheat(pcols,pver) ! heating rate due to snow melt (W/kg)
   real(r8), intent(out) :: precip(pcols)        ! rate of precipitation (kg / (m**2 * s))
   real(r8), intent(out) :: snowab(pcols)        ! rate of snow (kg / (m**2 * s))
   real(r8), intent(out) :: ice2pr(pcols,pver)   ! rate of conversion of ice to precip
   real(r8), intent(out) :: liq2pr(pcols,pver)   ! rate of conversion of liquid to precip
   real(r8), intent(out) :: liq2snow(pcols,pver) ! rate of conversion of liquid to snow
   real(r8), intent(out) :: rkflxprc(pcols,pverp)   ! grid-box mean RK flux_large_scale_cloud_rain+snow at interfaces (kg m^-2 s^-1)
   real(r8), intent(out) :: rkflxsnw(pcols,pverp)   ! grid-box mean RK flux_large_scale_cloud_snow at interfaces (kg m^-2 s^-1)
! intent(out)s here for pcond to pass to stratiform.F90 to be addflded/outflded
   real(r8), intent(out) :: pracwo(pcols,pver)      ! accretion of cloud water by rain (1/s)
   real(r8), intent(out) :: psacwo(pcols,pver)      ! accretion of cloud water by snow (1/s)
   real(r8), intent(out) :: psacio(pcols,pver)      ! accretion of cloud ice by snow (1/s)

   real(r8) nice2pr     ! rate of conversion of ice to snow
   real(r8) nliq2pr     ! rate of conversion of liquid to precip
   real(r8) nliq2snow   ! rate of conversion of liquid to snow
   real(r8) prodsnow(pcols,pver) ! rate of production of snow

!
! Local workspace
!
   real(r8) :: precab(pcols)        ! rate of precipitation (kg / (m**2 * s))
   integer i                 ! work variable
   integer iter              ! #iterations for precipitation calculation
   integer k                 ! work variable
   integer l                 ! work variable

   real(r8) cldm(pcols)          ! mean cloud fraction over the time step
   real(r8) cldmax(pcols)        ! max cloud fraction above
   real(r8) coef(pcols)          ! conversion time scale for condensate to rain
   real(r8) cwm(pcols)           ! cwat mixing ratio at midpoint of time step
   real(r8) cwn(pcols)           ! cwat mixing ratio at end
   real(r8) denom                ! work variable
   real(r8) dqsdt                ! change in sat spec. hum. wrt temperature
   real(r8) es(pcols)            ! sat. vapor pressure
   real(r8) fracw(pcols,pver)    ! relative importance of collection of liquid by rain
   real(r8) fsaci(pcols,pver)    ! relative importance of collection of ice by snow
   real(r8) fsacw(pcols,pver)    ! relative importance of collection of liquid by snow
   real(r8) fsaut(pcols,pver)    ! relative importance of ice auto conversion
   real(r8) fwaut(pcols,pver)    ! relative importance of warm cloud autoconversion
   real(r8) gamma(pcols)         ! d qs / dT
   real(r8) icwc(pcols)          ! in-cloud water content (kg/kg)
   real(r8) mincld               ! a small cloud fraction to avoid / zero
   real(r8) omeps                ! 1 minus epsilon
   real(r8),parameter ::omsm=0.99999_r8                 ! a number just less than unity (for rounding)
   real(r8) prprov(pcols)        ! provisional value of precip at btm of layer
   real(r8) prtmp                ! work variable
   real(r8) q(pcols,pver)        ! mixing ratio before time step ignoring condensate
   real(r8) qs(pcols)            ! spec. hum. of water vapor
   real(r8) qsn, esn             ! work variable
   real(r8) qsp(pcols,pver)      ! sat pt mixing ratio
   real(r8) qtl(pcols)           ! tendency which would saturate the grid box in deltat
   real(r8) qtmp, ttmp           ! work variable
   real(r8) relhum1(pcols)        ! relative humidity
   real(r8) relhum(pcols)        ! relative humidity
!!$   real(r8) tc                   ! crit temp of transition to ice
   real(r8) t(pcols,pver)        ! temp before time step ignoring condensate
   real(r8) tsp(pcols,pver)      ! sat pt temperature
   real(r8) pol                  ! work variable
   real(r8) cdt                  ! work variable
   real(r8) wtthick              ! work variable

! Extra local work space for cloud scheme modification       

   real(r8) cpohl                !Cp/Hlatv
   real(r8) hlocp                !Hlatv/Cp
   real(r8) dto2                 !0.5*deltat (delta=2.0*dt)
   real(r8) calpha(pcols)        !alpha of new C - E scheme formulation
   real(r8) cbeta (pcols)        !beta  of new C - E scheme formulation
   real(r8) cbetah(pcols)        !beta_hat at saturation portion 
   real(r8) cgamma(pcols)        !gamma of new C - E scheme formulation
   real(r8) cgamah(pcols)        !gamma_hat at saturation portion
   real(r8) rcgama(pcols)        !gamma/gamma_hat
   real(r8) csigma(pcols)        !sigma of new C - E scheme formulation
   real(r8) cmec1 (pcols)        !c1    of new C - E scheme formulation
   real(r8) cmec2 (pcols)        !c2    of new C - E scheme formulation
   real(r8) cmec3 (pcols)        !c3    of new C - E scheme formulation
   real(r8) cmec4 (pcols)        !c4    of new C - E scheme formulation
   real(r8) cmeres(pcols)        !residual cond of over-sat after cme and evapprec
   real(r8) ctmp                 !a scalar representation of cmeres
   real(r8) clrh2o               ! Ratio of latvap to water vapor gas const
   real(r8) ice(pcols,pver)    ! ice mixing ratio
   real(r8) liq(pcols,pver)    ! liquid mixing ratio
   real(r8) rcwn(pcols,2,pver), rliq(pcols,2,pver), rice(pcols,2,pver)
   real(r8) cwnsave(pcols,2,pver), cmesave(pcols,2,pver)
   real(r8) prodprecsave(pcols,2,pver)
   logical error_found
!
!------------------------------------------------------------
!
   clrh2o = hlatv/rh2o   ! Ratio of latvap to water vapor gas const
   omeps = 1.0_r8 - epsqs
#ifdef PERGRO
   mincld = 1.e-4_r8
   iter = 1   ! number of times to iterate the precipitation calculation
#else
   mincld = 1.e-4_r8
   iter = 2
#endif
!   omsm = 0.99999
   cpohl = cp/hlatv
   hlocp = hlatv/cp
   dto2=0.5_r8*deltat
!
! Constant for computing rate of evaporation of precipitation:
!
!!$   conke = 1.e-5
!!$   conke = 1.e-6
!
! initialize a few single level fields
!
   do i = 1,ncol
      precip(i) = 0.0_r8
      precab(i) = 0.0_r8
      snowab(i) = 0.0_r8
      cldmax(i) = 0.0_r8
   end do
!
! initialize multi-level fields 
!
   do k = 1,pver
      do i = 1,ncol
         q(i,k) = qn(i,k) 
         t(i,k) = tn(i,k)
!         q(i,k)=qn(i,k)-qtend(i,k)*deltat
!         t(i,k)=tn(i,k)-ttend(i,k)*deltat
    end do
   end do
   cme     (:ncol,:) = 0._r8
   evapprec(:ncol,:) = 0._r8
   prodprec(:ncol,:) = 0._r8
   evapsnow(:ncol,:) = 0._r8
   prodsnow(:ncol,:) = 0._r8
   evapheat(:ncol,:) = 0._r8
   meltheat(:ncol,:) = 0._r8
   prfzheat(:ncol,:) = 0._r8
   ice2pr(:ncol,:)   = 0._r8
   liq2pr(:ncol,:)   = 0._r8
   liq2snow(:ncol,:) = 0._r8
   fwaut(:ncol,:) = 0._r8
   fsaut(:ncol,:) = 0._r8
   fracw(:ncol,:) = 0._r8
   fsacw(:ncol,:) = 0._r8
   fsaci(:ncol,:) = 0._r8
   rkflxprc(:ncol,:) = 0._r8
   rkflxsnw(:ncol,:) = 0._r8

   pracwo(:ncol,:) = 0._r8
   psacwo(:ncol,:) = 0._r8
   psacio(:ncol,:) = 0._r8
!
! find the wet bulb temp and saturation value
! for the provisional t and q without condensation
!
   call findsp (lchnk, ncol, qn, tn, p, tsp, qsp)
   do 800 k = k1mb,pver
      call vqsatd (t(1,k), p(1,k), es, qs, gamma, ncol)
      do i = 1,ncol
         relhum(i) = q(i,k)/qs(i)
!
         cldm(i) = max(cldn(i,k),mincld)
!
! the max cloud fraction above this level
!
         cldmax(i) = max(cldmax(i), cldm(i))

! define the coefficients for C - E calculation

         calpha(i) = 1.0_r8/qs(i)
         cbeta (i) = q(i,k)/qs(i)**2*gamma(i)*cpohl
         cbetah(i) = 1.0_r8/qs(i)*gamma(i)*cpohl
         cgamma(i) = calpha(i)+hlatv*cbeta(i)/cp
         cgamah(i) = calpha(i)+hlatv*cbetah(i)/cp
         rcgama(i) = cgamma(i)/cgamah(i)

         if(cldm(i) > mincld) then
            icwc(i) = max(0._r8,cwat(i,k)/cldm(i))
         else
            icwc(i) = 0.0_r8
         endif
!PJR the above logic give zero icwc with nonzero cwat, dont like it!
!PJR generates problems with csigma
!PJR set the icwc to a very small number, so we can start from zero cloud cover and make some clouds
!         icwc(i) = max(1.e-8_r8,cwat(i,k)/cldm(i))

!
! initial guess of evaporation, will be updated within iteration
!
         evapprec(i,k) = conke*(1._r8 - cldm(i))*sqrt(precab(i)) &
                        *(1._r8 - min(relhum(i),1._r8))

!
! zero cmeres before iteration for each level
!
         cmeres(i)=0.0_r8

      end do
      do i = 1,ncol
!
! fractions of ice at this level
!
!!$         tc = t(i,k) - t0
!!$         fice(i,k) = max(0._r8,min(-tc*0.05,1.0_r8))
!
! calculate the cooling due to a phase change of the rainwater
! from above
!
         if (t(i,k) >= t0) then
            meltheat(i,k) =  -hlatf * snowab(i) * gravit/pdel(i,k)
            snowab(i) = 0._r8
         else
            meltheat(i,k) = 0._r8
         endif
      end do

!
! calculate cme and formation of precip. 
!
! The cloud microphysics is highly nonlinear and coupled with cme
! Both rain processes and cme are calculated iteratively.
! 
      do 100 l = 1,iter

        do i = 1,ncol

!
! calculation of cme has 4 scenarios
! ==================================
!
           if(relhum(i) > rhu00(i,k)) then
    
           ! 1. whole grid saturation
           ! ========================
              if(relhum(i) >= 0.999_r8 .or. cldm(i) >= 0.999_r8 ) then
                 cme(i,k)=(calpha(i)*qtend(i,k)-cbetah(i)*ttend(i,k))/cgamah(i)

           ! 2. fractional saturation
           ! ========================
              else
                 if (rhdfda(i,k) .eq. 0. .and. icwc(i).eq.0.) then
                    write (iulog,*) ' cldwat.F90:  empty rh cloud ', i, k, lchnk
                    write (iulog,*) ' relhum, iter ', relhum(i), l, rhu00(i,k), cldm(i), cldn(i,k)
                    call endrun ()
                 endif
                  csigma(i) = 1.0_r8/(rhdfda(i,k)+cgamma(i)*icwc(i))
                  cmec1(i) = (1.0_r8-cldm(i))*csigma(i)*rhdfda(i,k)
                  cmec2(i) = cldm(i)*calpha(i)/cgamah(i)+(1.0_r8-rcgama(i)*cldm(i))*   &
                             csigma(i)*calpha(i)*icwc(i)
                  cmec3(i) = cldm(i)*cbetah(i)/cgamah(i) +  &
                           (cbeta(i)-rcgama(i)*cldm(i)*cbetah(i))*csigma(i)*icwc(i)
                  cmec4(i) = csigma(i)*cgamma(i)*icwc(i)

                  ! Q=C-E=-C1*Al + C2*Aq - C3* At + C4*Er
  
                  cme(i,k) = -cmec1(i)*lctend(i,k) + cmec2(i)*qtend(i,k)  &
                             -cmec3(i)*ttend(i,k) + cmec4(i)*evapprec(i,k)
               endif

           ! 3. when rh < rhu00, evaporate existing cloud water
           ! ================================================== 
           else if(cwat(i,k) > 0.0_r8)then
              ! liquid water should be evaporated but not to exceed 
              ! saturation point. if qn > qsp, not to evaporate cwat
              cme(i,k)=-min(max(0._r8,qsp(i,k)-qn(i,k)),cwat(i,k))/deltat 

           ! 4. no condensation nor evaporation
           ! ==================================
           else
              cme(i,k)=0.0_r8
           endif

  
        end do    !end loop for cme update

! Because of the finite time step, 
! place a bound here not to exceed wet bulb point
! and not to evaporate more than available water
!
         do i = 1, ncol
            qtmp = qn(i,k) - cme(i,k)*deltat

! possibilities to have qtmp > qsp
!
!   1. if qn > qs(tn), it condenses; 
!      if after applying cme,  qtmp > qsp,  more condensation is applied. 
!      
!   2. if qn < qs, evaporation should not exceed qsp,
    
            if(qtmp > qsp(i,k)) then
              cme(i,k) = cme(i,k) + (qtmp-qsp(i,k))/deltat
            endif

!
! if net evaporation, it should not exceed available cwat
!
            if(cme(i,k) < -cwat(i,k)/deltat)  &
               cme(i,k) = -cwat(i,k)/deltat
!
! addition of residual condensation from previous step of iteration
!
            cme(i,k) = cme(i,k) + cmeres(i)

         end do

         !      limit cme for roundoff errors
         do i = 1, ncol
            cme(i,k) = cme(i,k)*omsm
         end do

         do i = 1,ncol
!
! as a safe limit, condensation should not reduce grid mean rh below rhu00
! 
           if(cme(i,k) > 0.0_r8 .and. relhum(i) > rhu00(i,k) )  &
              cme(i,k) = min(cme(i,k), (qn(i,k)-qs(i)*rhu00(i,k))/deltat)
!
! initial guess for cwm (mean cloud water over time step) if 1st iteration
!
           if(l < 2) then
             cwm(i) = max(cwat(i,k)+cme(i,k)*dto2,  0._r8)
           endif

         enddo

! provisional precipitation falling through model layer
         do i = 1,ncol
!!$            prprov(i) =  precab(i) + prodprec(i,k)*pdel(i,k)/gravit
! rain produced in this layer not too effective in collection process
            wtthick = max(0._r8,min(0.5_r8,((zi(i,k)-zi(i,k+1))/1000._r8)**2))
            prprov(i) =  precab(i) + wtthick*prodprec(i,k)*pdel(i,k)/gravit
         end do

! calculate conversion of condensate to precipitation by cloud microphysics 
         call findmcnew (lchnk   ,ncol    , &
                         k       ,prprov  ,snowab,  t       ,p        , &
                         cwm     ,cldm    ,cldmax  ,fice(1,k),coef    , &
                         fwaut(1,k),fsaut(1,k),fracw(1,k),fsacw(1,k),fsaci(1,k), &
                         seaicef, snowh, pracwo(1,k), psacwo(1,k), psacio(1,k))

!
! calculate the precip rate
!
         error_found = .false.
         do i = 1,ncol
            if (cldm(i) > 0) then  
!
! first predict the cloud water
!
               cdt = coef(i)*deltat
               if (cdt > 0.01_r8) then
                  pol = cme(i,k)/coef(i) ! production over loss
                  cwn(i) = max(0._r8,(cwat(i,k)-pol)*exp(-cdt)+ pol)
               else
                  cwn(i) = max(0._r8,(cwat(i,k) + cme(i,k)*deltat)/(1+cdt))
               endif
!
! now back out the tendency of net rain production
!
               prodprec(i,k) =  max(0._r8,cme(i,k)-(cwn(i)-cwat(i,k))/deltat)
            else
               prodprec(i,k) = 0.0_r8
               cwn(i) = 0._r8
            endif

            ! provisional calculation of conversion terms
            ice2pr(i,k) = prodprec(i,k)*(fsaut(i,k)+fsaci(i,k))
            liq2pr(i,k) = prodprec(i,k)*(fwaut(i,k)+fsacw(i,k)+fracw(i,k))
!old        liq2snow(i,k) = prodprec(i,k)*fsacw(i,k)

!           revision suggested by Jim McCaa
!           it controls the amount of snow hitting the sfc 
!           by forcing a lot of conversion of cloud liquid to snow phase
!           it might be better done later by an explicit representation of 
!           rain accreting ice (and freezing), or by an explicit freezing of raindrops
            liq2snow(i,k) = max(prodprec(i,k)*fsacw(i,k), fsnow(i,k)*liq2pr(i,k))

            ! bounds
            nice2pr = min(ice2pr(i,k),(cwat(i,k)+cme(i,k)*deltat)*fice(i,k)/deltat)
            nliq2pr = min(liq2pr(i,k),(cwat(i,k)+cme(i,k)*deltat)*(1._r8-fice(i,k))/deltat)
!            write(iulog,*) ' prodprec ', i, k, prodprec(i,k)
!            write(iulog,*) ' nliq2pr, nice2pr ', nliq2pr, nice2pr
            if (liq2pr(i,k).ne.0._r8) then
               nliq2snow = liq2snow(i,k)*nliq2pr/liq2pr(i,k)   ! correction
            else
               nliq2snow = liq2snow(i,k)
            endif

!           avoid roundoff problems generating negatives
            nliq2snow = nliq2snow*omsm
            nliq2pr = nliq2pr*omsm
            nice2pr = nice2pr*omsm
            
!           final estimates of conversion to precip and snow
            prodprec(i,k) = (nliq2pr + nice2pr)
            prodsnow(i,k) = (nice2pr + nliq2snow)

            rcwn(i,l,k) =  cwat(i,k) + (cme(i,k)-   prodprec(i,k))*deltat
            rliq(i,l,k) = (cwat(i,k) + cme(i,k)*deltat)*(1._r8-fice(i,k)) - nliq2pr * deltat
            rice(i,l,k) = (cwat(i,k) + cme(i,k)*deltat)* fice(i,k)      -    nice2pr                     *deltat

!           Save for sanity check later...  
!           Putting sanity checks inside loops 100 and 800 screws up the 
!           IBM compiler for reasons as yet unknown.  TBH
            cwnsave(i,l,k)      = cwn(i)
            cmesave(i,l,k)      = cme(i,k)
            prodprecsave(i,l,k) = prodprec(i,k)
!           End of save for sanity check later...  

!           final version of condensate to precip terms
            liq2pr(i,k) = nliq2pr
            liq2snow(i,k) = nliq2snow
            ice2pr(i,k) = nice2pr

            cwn(i) = rcwn(i,l,k)
!
! update any remaining  provisional values
!
            cwm(i) = (cwn(i) + cwat(i,k))*0.5_r8
!
! update in cloud water
!
            if(cldm(i) > mincld) then
               icwc(i) = cwm(i)/cldm(i)
            else
               icwc(i) = 0.0_r8
            endif
!PJR the above logic give zero icwc with nonzero cwat, dont like it!
!PJR generates problems with csigma
!PJR set the icwc to a very small number, so we can start from zero cloud cover and make some clouds
!         icwc(i) = max(1.e-8_r8,cwm(i)/cldm(i))

         end do              ! end of do i = 1,ncol

!
! calculate provisional value of cloud water for
! evaporation of precipitate (evapprec) calculation
!
      do i = 1,ncol
         qtmp = qn(i,k) - cme(i,k)*deltat
         ttmp = tn(i,k) + deltat/cp * ( meltheat(i,k)       &
              + (hlatv + hlatf*fice(i,k)) * cme(i,k) )
         esn = estblf(ttmp)
         qsn = min(epsqs*esn/(p(i,k) - omeps*esn),1._r8)
         qtl(i) = max((qsn - qtmp)/deltat,0._r8)
         relhum1(i) = qtmp/qsn
      end do
!
      do i = 1,ncol
#ifdef PERGRO
         evapprec(i,k) = conke*(1._r8 - max(cldm(i),mincld))* &
                         sqrt(precab(i))*(1._r8 - min(relhum1(i),1._r8))
#else
         evapprec(i,k) = conke*(1._r8 - cldm(i))*sqrt(precab(i)) &
                         *(1._r8 - min(relhum1(i),1._r8))
#endif
!
! limit the evaporation to the amount which is entering the box
! or saturates the box
!
         prtmp = precab(i)*gravit/pdel(i,k)
         evapprec(i,k) = min(evapprec(i,k), prtmp, qtl(i))*omsm
#ifdef PERGRO
!           zeroing needed for pert growth
         evapprec(i,k) = 0._r8
#endif
!
! Partition evaporation of precipitate between rain and snow using
! the fraction of snow falling into the box. Determine the heating
! due to evaporation. Note that evaporation is positive (loss of precip,
! gain of vapor) and that heating is negative.
         if (evapprec(i,k) > 0._r8) then
            evapsnow(i,k) = evapprec(i,k) * snowab(i) / precab(i)
            evapheat(i,k) = -hlatv * evapprec(i,k) - hlatf * evapsnow(i,k)
         else 
            evapsnow(i,k) = 0._r8
            evapheat(i,k) = 0._r8
         end if
! Account for the latent heat of fusion for liquid drops collected by falling snow
         prfzheat(i,k) = hlatf * liq2snow(i,k)
      end do

! now remove the residual of any over-saturation. Normally,
! the oversaturated water vapor should have been removed by 
! cme formulation plus constraints by wet bulb tsp/qsp
! as computed above. However, because of non-linearity,
! addition of (cme-evapprec) to update t and q may still cause
! a very small amount of over saturation. It is called a
! residual of over-saturation because theoretically, cme
! should have taken care of all of large scale condensation.
! 

       do i = 1,ncol
          qtmp = qn(i,k)-(cme(i,k)-evapprec(i,k))*deltat
          ttmp = tn(i,k) + deltat/cp * ( meltheat(i,k) + evapheat(i,k) + prfzheat(i,k)      &
              + (hlatv + hlatf*fice(i,k)) * cme(i,k) )
          esn = estblf(ttmp)
          qsn = min(epsqs*esn/(p(i,k) - omeps*esn),1._r8)
          !
          !Upper stratosphere and mesosphere, qsn calculated
          !above may be negative. Here just to skip it instead
          !of resetting it to 1 as in aqsat
          !
          if(qtmp > qsn .and. qsn > 0) then
             !calculate dqsdt, a more precise calculation
             !which taking into account different range of T 
             !can be found in aqsatd.F. Here follows
             !cond.F to calculate it.
             !
             denom = (p(i,k)-omeps*esn)*ttmp*ttmp
             dqsdt = clrh2o*qsn*p(i,k)/denom
             !
             !now extra condensation to bring air to just saturation
             !
             ctmp = (qtmp-qsn)/(1._r8+hlocp*dqsdt)/deltat
             cme(i,k) = cme(i,k)+ctmp
!
! save residual on cmeres to addtion to cme on entering next iteration
! cme exit here contain the residual but overrided if back to iteration
!
             cmeres(i) = ctmp
          else
             cmeres(i) = 0.0_r8
          endif
       end do
              
 100 continue              ! end of do l = 1,iter

!
! precipitation
!
      do i = 1,ncol
         precip(i) = precip(i) + pdel(i,k)/gravit * (prodprec(i,k) - evapprec(i,k))
         precab(i) = precab(i) + pdel(i,k)/gravit * (prodprec(i,k) - evapprec(i,k))
         if(precab(i).lt.0._r8) precab(i)=0._r8
!         snowab(i) = snowab(i) + pdel(i,k)/gravit * (prodprec(i,k)*fice(i,k) - evapsnow(i,k))
         snowab(i) = snowab(i) + pdel(i,k)/gravit * (prodsnow(i,k) - evapsnow(i,k))

         ! If temperature above freezing, all precip is rain flux.  if temperature below freezing, all precip is snow flux.
         rkflxprc(i,k+1) = precab(i)   !! making this consistent with other precip fluxes.  prc = rain + snow
         !!rkflxprc(i,k+1) = precab(i) - snowab(i)
         rkflxsnw(i,k+1) = snowab(i)

!!$         if ((precab(i)) < 1.e-10) then      
!!$            precab(i) = 0.
!!$            snowab(i) = 0.
!!$         endif
      end do
 800 continue                ! level loop (k=1,pver)

! begin sanity checks
   error_found = .false.
   do k = k1mb,pver
      do l = 1,iter
         do i = 1,ncol
	    if (abs(rcwn(i,l,k)).lt.1.e-300_r8) rcwn(i,l,k) = 0._r8
	    if (abs(rliq(i,l,k)).lt.1.e-300_r8) rliq(i,l,k) = 0._r8
	    if (abs(rice(i,l,k)).lt.1.e-300_r8) rice(i,l,k) = 0._r8
            if (rcwn(i,l,k).lt.0._r8) error_found = .true.
            if (rliq(i,l,k).lt.0._r8) error_found = .true.
            if (rice(i,l,k).lt.0._r8) error_found = .true.
         enddo
      enddo
   enddo
   if (error_found) then
      do k = k1mb,pver
         do l = 1,iter
            do i = 1,ncol
               if (rcwn(i,l,k).lt.0._r8) then
                  write(iulog,*) ' prob with neg rcwn1 ', rcwn(i,l,k),  &
                     cwnsave(i,l,k)
                  write(iulog,*) ' cwat, cme*deltat, prodprec*deltat ', &
                     cwat(i,k), cmesave(i,l,k)*deltat,               &
                     prodprecsave(i,l,k)*deltat,                     &
                     (cmesave(i,l,k)-prodprecsave(i,l,k))*deltat
                  call endrun('PCOND')
               endif
               if (rliq(i,l,k).lt.0._r8) then
                  write(iulog,*) ' prob with neg rliq1 ', rliq(i,l,k)
                  call endrun('PCOND')
               endif
               if (rice(i,l,k).lt.0._r8) then
                  write(iulog,*) ' prob with neg rice ', rice(i,l,k)
                  call endrun('PCOND')
               endif
            enddo
         enddo
      enddo
   end if
! end sanity checks

   return
end subroutine pcond

!##############################################################################

subroutine findmcnew (lchnk   ,ncol    , &
                      k       ,precab  ,snowab,  t       ,p       , &
                      cwm     ,cldm    ,cldmax  ,fice    ,coef    , &
                      fwaut   ,fsaut   ,fracw   ,fsacw   ,fsaci   , &
                      seaicef ,snowh,  pracwo, psacwo, psacio )
 
!----------------------------------------------------------------------- 
! 
! Purpose: 
! calculate the conversion of condensate to precipitate
! 
! Method: 
! See: Rasch, P. J, and J. E. Kristjansson, A Comparison of the CCM3
!  model climate using diagnosed and 
!  predicted condensate parameterizations, 1998, J. Clim., 11,
!  pp1587---1614.
! 
! Author: P. Rasch
! 
!-----------------------------------------------------------------------
   use phys_grid, only: get_rlat_all_p
   use comsrf,        only: landm
!
! input args
!
   integer, intent(in) :: lchnk                 ! chunk identifier
   integer, intent(in) :: ncol                  ! number of atmospheric columns
   integer, intent(in) :: k                     ! level index

   real(r8), intent(in) :: precab(pcols)        ! rate of precipitation from above (kg / (m**2 * s))
   real(r8), intent(in) :: t(pcols,pver)        ! temperature       (K)
   real(r8), intent(in) :: p(pcols,pver)        ! pressure          (Pa)
   real(r8), intent(in) :: cldm(pcols)          ! cloud fraction
   real(r8), intent(in) :: cldmax(pcols)        ! max cloud fraction above this level
   real(r8), intent(in) :: cwm(pcols)           ! condensate mixing ratio (kg/kg)
   real(r8), intent(in) :: fice(pcols)          ! fraction of cwat that is ice
   real(r8), intent(in) :: seaicef(pcols)       ! sea ice fraction 
   real(r8), intent(in) :: snowab(pcols)        ! rate of snow from above (kg / (m**2 * s))
    real(r8), intent(in) :: snowh(pcols)         ! Snow depth over land, water equivalent (m)

! output arguments
   real(r8), intent(out) :: coef(pcols)          ! conversion rate (1/s)
   real(r8), intent(out) :: fwaut(pcols)         ! relative importance of liquid autoconversion (a diagnostic)
   real(r8), intent(out) :: fsaut(pcols)         ! relative importance of ice autoconversion    (a diagnostic)
   real(r8), intent(out) :: fracw(pcols)         ! relative importance of rain accreting liquid (a diagnostic)
   real(r8), intent(out) :: fsacw(pcols)         ! relative importance of snow accreting liquid (a diagnostic)
   real(r8), intent(out) :: fsaci(pcols)         ! relative importance of snow accreting ice    (a diagnostic)
   real(r8), intent(out) :: pracwo(pcols)        ! accretion of cloud water by rain (1/s)
   real(r8), intent(out) :: psacwo(pcols)        ! accretion of cloud water by snow (1/s)
   real(r8), intent(out) :: psacio(pcols)        ! accretion of cloud ice by snow (1/s)


! work variables

   integer i
   integer ii
   integer ind(pcols)
   integer ncols

   real(r8), parameter :: degrad = 57.296_r8 ! divide by this to convert degrees to radians
   real(r8) capn                 ! local cloud particles / cm3
   real(r8) capnoice             ! local cloud particles when not over sea ice / cm3
   real(r8) ciaut                ! coefficient of autoconversion of ice (1/s)
   real(r8) cldloc(pcols)        ! non-zero amount of cloud
   real(r8) cldpr(pcols)         ! assumed cloudy volume occupied by rain and cloud
   real(r8) con1                 ! work constant
   real(r8) con2                 ! work constant
   real(r8) csacx                ! constant used for snow accreting liquid or ice
!!$   real(r8) dtice                ! interval for transition from liquid to ice
   real(r8) icemr(pcols)         ! in-cloud ice mixing ratio
   real(r8) icrit                ! threshold for autoconversion of ice
   real(r8) liqmr(pcols)         ! in-cloud liquid water mixing ratio
   real(r8) pracw                ! rate of rain accreting water
   real(r8) prlloc(pcols)        ! local rain flux in mm/day
   real(r8) prscgs(pcols)        ! local snow amount in cgs units
   real(r8) psaci                ! rate of collection of ice by snow (lin et al 1983)
   real(r8) psacw                ! rate of collection of liquid by snow (lin et al 1983)
   real(r8) psaut                ! rate of autoconversion of ice condensate
   real(r8) ptot                 ! total rate of conversion
   real(r8) pwaut                ! rate of autoconversion of liquid condensate
   real(r8) r3l                  ! volume radius
   real(r8) rainmr(pcols)        ! in-cloud rain mixing ratio
   real(r8) rat1                 ! work constant
   real(r8) rat2                 ! work constant
!!$   real(r8) rdtice               ! recipricol of dtice
   real(r8) rho(pcols)           ! density (mks units)
   real(r8) rhocgs               ! density (cgs units)
   real(r8) rlat(pcols)          ! latitude (radians)
   real(r8) snowfr               ! fraction of precipate existing as snow
   real(r8) totmr(pcols)         ! in-cloud total condensate mixing ratio
   real(r8) vfallw               ! fall speed of precipitate as liquid
   real(r8) wp                   ! weight factor used in calculating pressure dep of autoconversion
   real(r8) wsi                  ! weight factor for sea ice
   real(r8) wt                   ! fraction of ice
   real(r8) wland                ! fraction of land

!      real(r8) csaci
!      real(r8) csacw
!      real(r8) cwaut
!      real(r8) efact
!      real(r8) lamdas
!      real(r8) lcrit
!      real(r8) rcwm
!      real(r8) r3lc2
!      real(r8) snowmr(pcols)
!      real(r8) vfalls

   real(8) ftot

!     inline statement functions
   real(r8) heavy, heavym, a1, a2, heavyp, heavymp
   heavy(a1,a2) = max(0._r8,sign(1._r8,a1-a2))  ! heavyside function
   heavym(a1,a2) = max(0.01_r8,sign(1._r8,a1-a2))  ! modified heavyside function
!
! New heavyside functions to perhaps address error growth problems
!
   heavyp(a1,a2) = a1/(a2+a1+1.e-36_r8)
   heavymp(a1,a2) = (a1+0.01_r8*a2)/(a2+a1+1.e-36_r8)

!
! find all the points where we need to do the microphysics
! and set the output variables to zero
!
   ncols = 0
   do i = 1,ncol
      coef(i) = 0._r8
      fwaut(i) = 0._r8
      fsaut(i) = 0._r8
      fracw(i) = 0._r8
      fsacw(i) = 0._r8
      fsaci(i) = 0._r8
      liqmr(i) = 0._r8
      rainmr(i) = 0._r8
      if (cwm(i) > 1.e-20_r8) then
         ncols = ncols + 1
         ind(ncols) = i
      endif
   end do

!cdir nodep
!DIR$ CONCURRENT
   do ii = 1,ncols
      i = ind(ii)
!
! the local cloudiness at this level
!
      cldloc(i) = max(cldmin,cldm(i))
!
! a weighted mean between max cloudiness above, and this layer
!
      cldpr(i) = max(cldmin,(cldmax(i)+cldm(i))*0.5_r8)
!
! decompose the suspended condensate into
! an incloud liquid and ice phase component
!
      totmr(i) = cwm(i)/cldloc(i)
      icemr(i) = totmr(i)*fice(i)
      liqmr(i) = totmr(i)*(1._r8-fice(i))
!
! density
!
      rho(i) = p(i,k)/(287._r8*t(i,k))
      rhocgs = rho(i)*1.e-3_r8     ! density in cgs units
!
! decompose the precipitate into a liquid and ice phase
!
      if (t(i,k) > t0) then
         vfallw = convfw/sqrt(rho(i))
         rainmr(i) = precab(i)/(rho(i)*vfallw*cldpr(i))
         snowfr = 0
!        snowmr(i)
      else
         snowfr = 1
         rainmr(i) = 0._r8
      endif
!     rainmr(i) = (precab(i)-snowab(i))/(rho(i)*vfallw*cldpr(i))
!
! local snow amount in cgs units
!
      prscgs(i) = precab(i)/cldpr(i)*0.1_r8*snowfr
!     prscgs(i) = snowab(i)/cldpr(i)*0.1
!
! local rain amount in mm/day
!
      prlloc(i) = precab(i)*86400._r8/cldpr(i)
   end do

   con1 = 1._r8/(1.333_r8*pi)**0.333_r8 * 0.01_r8 ! meters
!
! calculate the conversion terms
!
   call get_rlat_all_p(lchnk, ncol, rlat)

!cdir nodep
!DIR$ CONCURRENT
   do ii = 1,ncols
      i = ind(ii)
      rhocgs = rho(i)*1.e-3_r8     ! density in cgs units
!
! exponential temperature factor
!
!        efact = exp(0.025*(t(i,k)-t0))
!
! some temperature dependent constants
!
!!$      wt = min(1._r8,max(0._r8,(t0-t(i,k))*rdtice))
      wt = fice(i)
      icrit = icritc*wt + icritw*(1-wt)
!
! jrm Reworked droplet number concentration algorithm
      ! Start with pressure-dependent value appropriate for continental air
      ! Note: reltab has a temperature dependence here
      capn = capnw + (capnc-capnw) * min(1._r8,max(0._r8,1.0_r8-(p(i,k)-0.8_r8*p(i,pver))/(0.2_r8*p(i,pver))))
      ! Modify for snow depth over land
      capn = capn + (capnc-capn) * min(1.0_r8,max(0.0_r8,snowh(i)*10._r8))
      ! Ramp between polluted value over land to clean value over ocean.
      capn = capn + (capnc-capn) * min(1.0_r8,max(0.0_r8,1.0_r8-landm(i,lchnk)))
      ! Ramp between the resultant value and a sea ice value in the presence of ice.
      capn = capn + (capnsi-capn) * min(1.0_r8,max(0.0_r8,seaicef(i)))
! end jrm
!      
#ifdef DEBUG2
      if ( (lat(i) == latlook(1)) .or. (lat(i) == latlook(2)) ) then
         if (i == ilook(1)) then
            write(iulog,*) ' findmcnew: lat, k, seaicef, landm, wp, capnoice, capn ', &
                 lat(i), k, seaicef(i), landm(i,lat(i)), wp, capnoice, capn
         endif
      endif
#endif

!
! useful terms in following calculations
!
      rat1 = rhocgs/rhow
      rat2 = liqmr(i)/capn
      con2 = (rat1*rat2)**0.333_r8
!
! volume radius
!
!        r3l = (rhocgs*liqmr(i)/(1.333*pi*capn*rhow))**0.333 * 0.01 ! meters
      r3l = con1*con2
!
! critical threshold for autoconversion if modified for mixed phase
! clouds to mimic a bergeron findeisen process
! r3lc2 = r3lcrit*(1.-0.5*fice(i)*(1-fice(i)))
!
! autoconversion of liquid
!
!        cwaut = 2.e-4
!        cwaut = 1.e-3
!        lcrit = 2.e-4
!        lcrit = 5.e-4
!        pwaut = max(0._r8,liqmr(i)-lcrit)*cwaut
!
! pwaut is following tripoli and cotton (and many others)
! we reduce the autoconversion below critpr, because these are regions where
! the drop size distribution is likely to imply much smaller collector drops than
! those relevant for a cloud distribution corresponding to the value of effc = 0.55
! suggested by cotton (see austin 1995 JAS, baker 1993)

! easy to follow form
!        pwaut = capc*liqmr(i)**2*rhocgs/rhow
!    $           *(liqmr(i)*rhocgs/(rhow*capn))**(.333)
!    $           *heavy(r3l,r3lcrit)
!    $           *max(0.10_r8,min(1._r8,prlloc(i)/critpr))
! somewhat faster form
#define HEAVYNEW
#ifdef HEAVYNEW
!#ifdef PERGRO
      pwaut = capc*liqmr(i)**2*rat1*con2*heavymp(r3l,r3lcrit) * &
              max(0.10_r8,min(1._r8,prlloc(i)/critpr))
#else
      pwaut = capc*liqmr(i)**2*rat1*con2*heavym(r3l,r3lcrit)* &
              max(0.10_r8,min(1._r8,prlloc(i)/critpr))
#endif
!
! autoconversion of ice
!
!        ciaut = ciautb*efact
      ciaut = ciautb
!        psaut = capc*totmr(i)**2*rhocgs/rhoi
!     $           *(totmr(i)*rhocgs/(rhoi*capn))**(.333)
!
! autoconversion of ice condensate
!
#ifdef PERGRO
      psaut = heavyp(icemr(i),icrit)*icemr(i)*ciaut
#else
      psaut = max(0._r8,icemr(i)-icrit)*ciaut
#endif
!
! collection of liquid by rain
!
!        pracw = cracw*rho(i)*liqmr(i)*rainmr(i) !(beheng 1994)
      pracw = cracw*rho(i)*sqrt(rho(i))*liqmr(i)*rainmr(i) !(tripoli and cotton)

      pracwo(i)=pracw

!!      pracw = 0.
!
! the following lines calculate the slope parameter and snow mixing ratio
! from the precip rate using the equations found in lin et al 83
! in the most natural form, but it is expensive, so after some tedious
! algebraic manipulation you can use the cheaper form found below
!            vfalls = c*gam4pd/(6*lamdas**d)*sqrt(rhonot/rhocgs)
!     $               *0.01   ! convert from cm/s to m/s
!            snowmr(i) = snowfr*precab(i)/(rho(i)*vfalls*cldpr(i))
!            snowmr(i) = ( prscgs(i)*mcon02 * (rhocgs**mcon03) )**mcon04
!            lamdas = (prhonos/max(rhocgs*snowmr(i),small))**0.25
!            csacw = mcon01*sqrt(rhonot/rhocgs)/(lamdas**thrpd)
!
! coefficient for collection by snow independent of phase
!
      csacx = mcon07*rhocgs**mcon08*prscgs(i)**mcon05

!
! collection of liquid by snow (lin et al 1983)
!
      psacw = csacx*liqmr(i)*esw
#ifdef PERGRO
! this is necessary for pergro
      psacw = 0._r8
#endif

      psacwo(i)=psacw

!
! collection of ice by snow (lin et al 1983)
!
      psaci = csacx*icemr(i)*esi
!
      psacio(i)=psaci

! total conversion of condensate to precipitate
!
      ptot = pwaut + psaut + pracw + psacw + psaci
!
! the recipricol of cloud water amnt (or zero if no cloud water)
!
!         rcwm =  totmr(i)/(max(totmr(i),small)**2)
!
! turn the tendency back into a loss rate (1/seconds)
!
      if (totmr(i) > 0._r8) then
         coef(i) = ptot/totmr(i)
      else
         coef(i) = 0._r8
      endif

      if (ptot.gt.0._r8) then
         fwaut(i) = pwaut/ptot
         fsaut(i) = psaut/ptot
         fracw(i) = pracw/ptot
         fsacw(i) = psacw/ptot
         fsaci(i) = psaci/ptot
      else
         fwaut(i) = 0._r8
         fsaut(i) = 0._r8
         fracw(i) = 0._r8
         fsacw(i) = 0._r8
         fsaci(i) = 0._r8
      endif

      ftot = fwaut(i)+fsaut(i)+fracw(i)+fsacw(i)+fsaci(i)
!      if (abs(ftot-1._r8).gt.1.e-14_r8.and.ftot.ne.0._r8) then
!         write(iulog,*) ' something is wrong in findmcnew ', ftot, &
!              fwaut(i),fsaut(i),fracw(i),fsacw(i),fsaci(i)
!         write(iulog,*) ' unscaled ', ptot, &
!              pwaut,psaut,pracw,psacw,psaci
!         write(iulog,*) ' totmr, liqmr, icemr ', totmr(i), liqmr(i), icemr(i)
!         call endrun()
!      endif
   end do
#ifdef DEBUG
   i = icollook(nlook)
   if (lchnk == lchnklook(nlook) ) then
      write(iulog,*)
      write(iulog,*) '------', k, i, lchnk
      write(iulog,*) ' liqmr, rainmr,precab ', liqmr(i), rainmr(i), precab(i)*8.64e4_r8
      write(iulog,*) ' frac: waut,saut,racw,sacw,saci ', &
           fwaut(i), fsaut(i), fracw(i), fsacw(i), fsaci(i)
   endif
#endif

   return
end subroutine findmcnew

!##############################################################################

subroutine findsp (lchnk, ncol, q, t, p, tsp, qsp)
!----------------------------------------------------------------------- 
! 
! Purpose: 
!     find the wet bulb temperature for a given t and q
!     in a longitude height section
!     wet bulb temp is the temperature and spec humidity that is 
!     just saturated and has the same enthalpy
!     if q > qs(t) then tsp > t and qsp = qs(tsp) < q
!     if q < qs(t) then tsp < t and qsp = qs(tsp) > q
!
! Method: 
! a Newton method is used
! first guess uses an algorithm provided by John Petch from the UKMO
! we exclude points where the physical situation is unrealistic
! e.g. where the temperature is outside the range of validity for the
!      saturation vapor pressure, or where the water vapor pressure
!      exceeds the ambient pressure, or the saturation specific humidity is 
!      unrealistic
! 
! Author: P. Rasch
! 
!-----------------------------------------------------------------------
!
!     input arguments
!
   integer, intent(in) :: lchnk                 ! chunk identifier
   integer, intent(in) :: ncol                  ! number of atmospheric columns

   real(r8), intent(in) :: q(pcols,pver)        ! water vapor (kg/kg)
   real(r8), intent(in) :: t(pcols,pver)        ! temperature (K)
   real(r8), intent(in) :: p(pcols,pver)        ! pressure    (Pa)
!
! output arguments
!
   real(r8), intent(out) :: tsp(pcols,pver)      ! saturation temp (K)
   real(r8), intent(out) :: qsp(pcols,pver)      ! saturation mixing ratio (kg/kg)
!
! local variables
!
   integer i                 ! work variable
   integer k                 ! work variable
   logical lflg              ! work variable
   integer iter              ! work variable
   integer l                 ! work variable
   logical :: error_found

   real(r8) omeps                ! 1 minus epsilon
   real(r8) trinv                ! work variable
   real(r8) es                   ! sat. vapor pressure
   real(r8) desdt                ! change in sat vap pressure wrt temperature
!     real(r8) desdp                ! change in sat vap pressure wrt pressure
   real(r8) dqsdt                ! change in sat spec. hum. wrt temperature
   real(r8) dgdt                 ! work variable
   real(r8) g                    ! work variable
   real(r8) weight(pcols)        ! work variable
   real(r8) hlatsb               ! (sublimation)
   real(r8) hlatvp               ! (vaporization)
   real(r8) hltalt(pcols,pver)   ! lat. heat. of vap.
   real(r8) tterm                ! work var.
   real(r8) qs                   ! spec. hum. of water vapor
   real(r8) tc                   ! crit temp of transition to ice

! work variables
   real(r8) t1, q1, dt, dq
   real(r8) dtm, dqm
   real(r8) qvd, a1, tmp
   real(r8) rair
   real(r8) r1b, c1, c2, c3
   real(r8) denom
   real(r8) dttol
   real(r8) dqtol
   integer doit(pcols) 
   real(r8) enin(pcols), enout(pcols)
   real(r8) tlim(pcols)

   omeps = 1.0_r8 - epsqs
   trinv = 1.0_r8/ttrice
   a1 = 7.5_r8*log(10._r8)
   rair =  287.04_r8
   c3 = rair*a1/cp
   dtm = 0._r8    ! needed for iter=0 blowup with f90 -ei
   dqm = 0._r8    ! needed for iter=0 blowup with f90 -ei
   dttol = 1.e-4_r8 ! the relative temp error tolerance required to quit the iteration
   dqtol = 1.e-4_r8 ! the relative moisture error tolerance required to quit the iteration
!  tmin = 173.16 ! the coldest temperature we can deal with
!
! max number of times to iterate the calculation
   iter = 8
!
   do k = k1mb,pver

!
! first guess on the wet bulb temperature
!
      do i = 1,ncol

#ifdef DEBUG
         if ( (lchnk == lchnklook(nlook) ) .and. (i == icollook(nlook) ) ) then
            write(iulog,*) ' '
            write(iulog,*) ' level, t, q, p', k, t(i,k), q(i,k), p(i,k)
         endif
#endif
! limit the temperature range to that relevant to the sat vap pres tables
#if ( ! defined WACCM_PHYS )
         tlim(i) = min(max(t(i,k),173._r8),373._r8)
#else
         tlim(i) = min(max(t(i,k),128._r8),373._r8)
#endif
         es = estblf(tlim(i))
         denom = p(i,k) - omeps*es
         qs = epsqs*es/denom
         doit(i) = 0
         enout(i) = 1._r8
! make sure a meaningful calculation is possible
         if (p(i,k) > 5._r8*es .and. qs > 0._r8 .and. qs < 0.5_r8) then
!
! Saturation specific humidity
!
             qs = min(epsqs*es/denom,1._r8)
!
! "generalized" analytic expression for t derivative of es
!  accurate to within 1 percent for 173.16 < t < 373.16
!
! Weighting of hlat accounts for transition from water to ice
! polynomial expression approximates difference between es over
! water and es over ice from 0 to -ttrice (C) (min of ttrice is
! -40): required for accurate estimate of es derivative in transition
! range from ice to water also accounting for change of hlatv with t
! above freezing where const slope is given by -2369 j/(kg c) = cpv - cw
!
             tc     = tlim(i) - t0
             lflg   = (tc >= -ttrice .and. tc < 0.0_r8)
             weight(i) = min(-tc*trinv,1.0_r8)
             hlatsb = hlatv + weight(i)*hlatf
             hlatvp = hlatv - 2369.0_r8*tc
             if (tlim(i) < t0) then
                hltalt(i,k) = hlatsb
             else
                hltalt(i,k) = hlatvp
             end if
             enin(i) = cp*tlim(i) + hltalt(i,k)*q(i,k)

! make a guess at the wet bulb temp using a UKMO algorithm (from J. Petch)
             tmp =  q(i,k) - qs
             c1 = hltalt(i,k)*c3
             c2 = (tlim(i) + 36._r8)**2
             r1b    = c2/(c2 + c1*qs)
             qvd   = r1b*tmp
             tsp(i,k) = tlim(i) + ((hltalt(i,k)/cp)*qvd)
#ifdef DEBUG
             if ( (lchnk == lchnklook(nlook) ) .and. (i == icollook(nlook) ) ) then
                write(iulog,*) ' relative humidity ', q(i,k)/qs
                write(iulog,*) ' first guess ', tsp(i,k)
             endif
#endif
             es = estblf(tsp(i,k))
             qsp(i,k) = min(epsqs*es/(p(i,k) - omeps*es),1._r8)
          else
             doit(i) = 1
             tsp(i,k) = tlim(i)
             qsp(i,k) = q(i,k)
             enin(i) = 1._r8
          endif
       end do   ! end do i
!
! now iterate on first guess
!
      do l = 1, iter
         dtm = 0
         dqm = 0
         do i = 1,ncol
            if (doit(i) == 0) then
               es = estblf(tsp(i,k))
!
! Saturation specific humidity
!
               qs = min(epsqs*es/(p(i,k) - omeps*es),1._r8)
!
! "generalized" analytic expression for t derivative of es
! accurate to within 1 percent for 173.16 < t < 373.16
!
! Weighting of hlat accounts for transition from water to ice
! polynomial expression approximates difference between es over
! water and es over ice from 0 to -ttrice (C) (min of ttrice is
! -40): required for accurate estimate of es derivative in transition
! range from ice to water also accounting for change of hlatv with t
! above freezing where const slope is given by -2369 j/(kg c) = cpv - cw
!
               tc     = tsp(i,k) - t0
               lflg   = (tc >= -ttrice .and. tc < 0.0_r8)
               weight(i) = min(-tc*trinv,1.0_r8)
               hlatsb = hlatv + weight(i)*hlatf
               hlatvp = hlatv - 2369.0_r8*tc
               if (tsp(i,k) < t0) then
                  hltalt(i,k) = hlatsb
               else
                  hltalt(i,k) = hlatvp
               end if
               if (lflg) then
                  tterm = pcf(1) + tc*(pcf(2) + tc*(pcf(3)+tc*(pcf(4) + tc*pcf(5))))
               else
                  tterm = 0.0_r8
               end if
               desdt = hltalt(i,k)*es/(rgasv*tsp(i,k)*tsp(i,k)) + tterm*trinv
               dqsdt = (epsqs + omeps*qs)/(p(i,k) - omeps*es)*desdt
!              g = cp*(tlim(i)-tsp(i,k)) + hltalt(i,k)*q(i,k)- hltalt(i,k)*qsp(i,k)
               g = enin(i) - (cp*tsp(i,k) + hltalt(i,k)*qsp(i,k))
               dgdt = -(cp + hltalt(i,k)*dqsdt)
               t1 = tsp(i,k) - g/dgdt
               dt = abs(t1 - tsp(i,k))/t1
               tsp(i,k) = max(t1,tmin)
               es = estblf(tsp(i,k))
               q1 = min(epsqs*es/(p(i,k) - omeps*es),1._r8)
               dq = abs(q1 - qsp(i,k))/max(q1,1.e-12_r8)
               qsp(i,k) = q1
#ifdef DEBUG
               if ( (lchnk == lchnklook(nlook) ) .and. (i == icollook(nlook) ) ) then
                  write(iulog,*) ' rel chg lev, iter, t, q ', k, l, dt, dq, g
               endif
#endif
               dtm = max(dtm,dt)
               dqm = max(dqm,dq)
! if converged at this point, exclude it from more iterations
               if (dt < dttol .and. dq < dqtol) then
                  doit(i) = 2
               endif
               enout(i) = cp*tsp(i,k) + hltalt(i,k)*qsp(i,k)
! bail out if we are too near the end of temp range
#if ( ! defined WACCM_PHYS )
               if (tsp(i,k) < 174.16_r8) then
#else
               if (tsp(i,k) < 130.16_r8) then
#endif
                  doit(i) = 4
               endif
            else
            endif
         end do              ! do i = 1,ncol

         if (dtm < dttol .and. dqm < dqtol) then
            go to 10
         endif

      end do                 ! do l = 1,iter
10    continue

      error_found = .false.
      if (dtm > dttol .or. dqm > dqtol) then
         do i = 1,ncol
            if (doit(i) == 0) error_found = .true.
         end do
         if (error_found) then
            do i = 1,ncol
               if (doit(i) == 0) then
                  write(iulog,*) ' findsp not converging at point i, k ', i, k
                  write(iulog,*) ' t, q, p, enin ', t(i,k), q(i,k), p(i,k), enin(i)
                  write(iulog,*) ' tsp, qsp, enout ', tsp(i,k), qsp(i,k), enout(i)
                  call endrun ('FINDSP')
               endif
            end do
         endif
      endif
      do i = 1,ncol
         if (doit(i) == 2 .and. abs((enin(i)-enout(i))/(enin(i)+enout(i))) > 1.e-4_r8) then
            error_found = .true.
         endif
      end do
      if (error_found) then
         do i = 1,ncol
            if (doit(i) == 2 .and. abs((enin(i)-enout(i))/(enin(i)+enout(i))) > 1.e-4_r8) then
               write(iulog,*) ' the enthalpy is not conserved for point ', &
                  i, k, enin(i), enout(i)
               write(iulog,*) ' t, q, p, enin ', t(i,k), q(i,k), p(i,k), enin(i)
               write(iulog,*) ' tsp, qsp, enout ', tsp(i,k), qsp(i,k), enout(i)
               call endrun ('FINDSP')
            endif
         end do
      endif
      
   end do                    ! level loop (k=1,pver)

   return
end subroutine findsp

end module cldwat
